On the Mechanical Interaction of Light With Homogeneous Liquids 



o 

(N 
X 



Michiel de Reus and Neil V. Budko 

Numerical Analysis, Delft Institute of Applied Mathematics, 
Faculty of Electrical Engineering, Mathematics and Computer Science, 
Delft University of Technology, Mekelweg 4, 2628 CD Delft, The Netherlands 

(Dated: May 2, 2013) 

We investigate one of the consequences of the three competing models describing the mechanical 
interaction of light with a dielectric medium. According to both the Abraham and Minkowski models 
the time-averaged force density is zero inside a homogeneous dielectric, whereas the induced-current 
Lorentz force model predicts a non-zero force density. We argue that the latter force, if exists, could 
drive a hydrodynamic flow inside a homogeneous fluid. Our numerical experiments show that such 
flows have distinct spatial patterns and may influence the dynamics of particles in a water-based 
single-beam optical trap. 
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I. INTRODUCTION 

The mechanical force exerted by light on dielectric ob- 
jects enables capturing, manipulation, and analysis of 
particles at the micro- and nano-meter scales 0, Q . Of- 
ten, especially in biological applications, the experimen- 
tal set-up of an optical trap involves a liquid host medium 
Q . The dynamics of particles in such traps is influenced 
not only by the optical and Brownian forces but by the 
drag force of the liquid as well. If operating in the neigh- 
borhood of an absorption band, a significant thermal flow 
may be induced by the light itself, and needs to be taken 
into account [4|-|(|. The analysis of the host- medium 
hydrodynamics is important for understanding the ob- 
served non-conservative motion of trapped particles Q, 
which at the moment is mostly attributed to the pres- 
ence of the so-called scattering force, Brownian 'motors' 
[H, and thermal/thermophoretic flows. The emergent 
field of optofluidics also requires a better understanding 
of the mechanical interaction of light with liquids [ij-[l2j . 

A question arises whether absorption/extinction is the 
only mechanism by which a light-driven hydrodynamic 
flow may be induced. Can a flow be created in a loss- 
less liquid? Since the host liquid is typically optically 
homogeneous the question can be narrowed down to the 
mechanical forces exerted by the electromagnetic field in- 
side a homogeneous dielectric. Yet, no simple answer is 
available at the moment. The two dominant approaches, 
associated with the names of Abraham and Minkowski, 
although differing on the form of the mechanical momen- 
tum of the electromagnetic field, seem to agree that the 
time-averaged force density (the Helmholtz force density) 
in a homogeneous medium should be exactly zero [13| . 
On the contrary, the third approach, that treats the in- 
duced (polarization) currents on equal footing with the 
externally imposed ones, leads to a non-zero Lorentz force 
density. Thus, one faces a fundamental question: does 
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the electromagnetic field exert a mechanical force on the 
polarization currents induced by that same field? 

It has been shown that all three approaches give the 
same total time-averaged mechanical force on a finite ho- 
mogeneous solid body, even though the Helmholtz force 
density contributes only at the boundary surface while 
the Lorentz force density acts across the whole volume 
of the object [2]. However, recently it has been sug- 
gested that a deformable homogeneous body would re- 
spond differently to these forces 13]. Although, the in- 
ternal stresses induced by light in thin films may appear 
to be difficult to detect experimentally, if such internal 
stresses do exist, they might induce hydrodynamic flows 
that are much easier to observe and should, in fact, be 
taken into account during the optical trapping. 

The quantitative analysis of light-driven flows for a 
typical single-beam optical trap presented here is aimed 
at working out the consequences of the induced-current 
Lorentz force density model. We consider both the trans- 
parent and absorptive spectral regions. Our simulations 
demonstrate that according to this model the light-driven 
flows should be significant (in the order of tens of microm- 
eters per second) even with small laser powers (~ 1 mW). 
Such flows have very distinctive patterns that could be 
detected in an experiment. 

The paper is organized as follows. First, we briefly de- 
scribe the three forms of the electromagnetic momentum 
conservation law and the associated time-averaged force 
densities. We devote a separate section to the question 
of modeling of the electromagnetic field in a single-beam 
optical trap, where we show that a standard Gaussian 
beam model is not applicable for the purposes of force 
computations as it fails to satisfy the momentum conser- 
vation law. We resort to an alternative method that mim- 
ics the Gaussian beam by a linear superposition of exact 
fundamental solutions of the Maxwell equations. Then, 
we introduce the mechanical coupling between the elec- 
tromagnetic field and the host liquid via the body force 
term of the incompressible Navier-Stokes equations and 
the Boussinesq approximation for the thermally-driven 
flows. Finally, we present the results of numerical exper- 
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iments and our conclusions. 



II. MODELS OF THE ELECTROMAGNETIC 
FORCE DENSITY 

There exist two mathematically equivalent views on 
the electromagnetic field in media. In vacuum the 
Maxwell equations have the form: 



-V xH + eoftE = -J oxt , 
V x E + nodtH = -Q, 



(1) 



with E(x, t), H(x, t) denoting the electric and magnetic 
field strengths, and J oxt (x, t) - the external (i.e., inde- 
pendent of the field) electric current density. 

A medium different from vacuum can be introduced 
either via the concept of electric and magnetic fluxes D 
and B entering the Maxwell equations as 



-V x H + d t T> = J° 
V x E + d t B = -0, 



(2) 



or via the induced electric and magnetic current densities 
J ind and K ind as 



-V x H + sM = -J cxt - J' 
V x E + /! <9 t H = -K ind , 



(3) 



Supplied with appropriate constitutive relations, for ex- 
ample, 



D = eE, B = fiH, 



(e - e )5 t E, K ind = (/U — Ho)d t H, 



(4) 



both formulations lead to exactly the same Maxwell's 
equations for the fields E and H. 

In general, the momentum conservation law has the 
form: 



v • t - ap = f , 



(5) 



where T is the stress tensor density, P is the momen- 
tum density, and f is the force density, which is our main 
concern. The actual form of this law depends on the 
definition and interpretation of the terms. For example, 
in the lossless case described by the constitutive relations 
(HI the Abraham and Minkowski expressions for the force 
density in the part of the domain without external cur- 
rents/charges are 



where 



f A = f H + ^V^9 t (ExH), 

Cr 

f M _ f H 



f H = -i(E.E)V£-i(H-H)V/i 



(6) 
(7) 

(8) 



is the Helmholtz force density. In a general medium this 
force density is defined as: 

f H = i (D • VE - E • VD) + i (B • VH - H • VB) . 

(9) 

Obviously, after time averaging over a period of harmonic 
oscillations the Helmholtz force will be the only non-zero 
contribution to both the Abraham and Minkowski force 
densities and it will be zero in a homogeneous lossless 
medium, i.e., 

(f A ) = (f M ) = (f H ) = 0, e,n= const. (10) 

The Abraham and Minkowski expressions follow from the 
fluxes-based approach to the medium ([2]). If, on the other 
hand, one starts with the induced-currents approach , 
then the force density turns out to be: 

^ind 



Pi 



{ L = p md E + p ma H + ^JHm x jj £oK md x j, 
t pt 

V-J ind dt', p™ d = / V-K ind dt' 



ind 
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In a lossless homogeneous medium this generalized 
Lorentz force density reduces to: 

f L = n (e - e ){d t E) x H - e (M - Mo)(d t H) x E. 

(12) 

Assuming time-harmonic fields with the angular fre- 
quency uj and a linear non-magnetic homogeneous dis- 
persive medium we perform the averaging of the general 
expression (|11[) over a period of oscillations and arrive at 
the following result: 
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<f L ) = -wpo(e' - £o)Im{S} + -w/z e"Re{S}, (13) 
where the complex Poynting vector is defined as 

S = ExH*, (14) 

Here, the complex field amplitudes E, H satisfy the 
frequency-domain Maxwell's equations, (...)* denotes 
the complex conjugation, and e(uj) = e'(uj) + ie"(uj) is 
the complex permittivity of the host liquid. 

The heat produced by the electromagnetic field is given 
by the following standard time-averaged expression de- 
rived from the energy conservation law: 



q = -u;e"E ■ E* 
2 



III. SINGLE-BEAM OPTICAL TRAP 



(15) 



In optics it is common to model the spatial distribu- 
tion of the field of a single-beam optical trap as a tightly 
focused Gaussian beam. Unfortunately, this model is 
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FIG. 1: Normalized residuals of the Maxwell equations (left) and the momentum conservation law (right) as a function of the 
beam parameter \s\, see eq. (|16p. for different discretization levels (n = 32,64,96, 128 uniformly spaced grid-points along each 
of the 10 /im sides of a cube). Solid lines - dipole-based model, dashed lines - Gaussian beam approximation. As the beam 
waist gets wider (smaller |s| values) the residual of the Maxwell equations for the Gaussian approximation becomes smaller and 
is limited by the discretization error. The residual of the momentum conservation law stays relatively high irrespectively of the 
discretization level for the Gaussian approximation, whereas it is consistently smaller and improves with finer discretization 
levels for the dipole model of the beam. 



not really suitable for the calculation of the force density 
since it violates the conservation of momentum to an ex- 
tent that may invalidate the results of numerical experi- 
ments. Indeed, the expression (fTBf represents the electro- 
magnetic force density only if the field quantities appear- 
ing in it are the exact solutions of the frequency-domain 
versions of the Maxwell equations ([2]) or ([3]) , whereas the 
Gaussian beam is not an exact solution of the Maxwell 
equations. 

Figure [1] (left) gives the norm of the residual of the 
frequency-domain Maxwell's equations computed numer- 
ically via the finite-volume approximation of the equa- 
tions on different progressively refined grids. The hori- 
zontal axis is the Gaussian beam parameter 



where Ao is the vacuum wavelength of light, n is the com- 
plex refractive index of the medium, and wq is the beam 
waist. Thus, the larger is the relative beam waist wq/Xo, 
the smaller is \s\. The residual of the Gaussian approxi- 
mation (dashed lines) reduces for larger waists, generally 
following the C(|s| 2 ) trend - the order of the approxima- 
tion. It also improves with the grid refinement. This, 
however, is merely a consequence of the finite- volume 
approximation of the spatial derivatives in the Maxwell 
equations. 

Figure [1] (right) shows the residual of the frequency- 
domain version of the momentum conservation law (|5|), 
where the left-hand side was estimated using the numer- 
ical surface integration over the elementary cells. Here 
too the Gaussian approximation (dashed lines) improves 
with the increase in the beam waist. However, refining 
the discretization does not help any more. Hence, we 



conclude that the error is mainly due to the analytical 
mismatch in the spatial structure of the fields. 

A way to improve the single-beam optical trap model 
is to represent the electromagnetic field as a superposi- 
tion of exact analytical solutions of the Maxwell equa- 
tions due to elementary dipole sources. With the help of 
these fundamental solutions we can model a focused laser 
beam by distributing an array of dipoles over a plane just 
outside the computational domain and tuning their am- 
plitudes, phases, and polarizations (dipole moments) so 
that they reproduce the spatial pattern of the electric 
field of an ideal focused Gaussian beam passing through 
that plane. Basically, this can be viewed as a numerical 
implementation of the Huygens principle. 

Solid lines in Figure [T] demonstrate that the residuals 
for the dipole-based model of the beam stay roughly con- 
stant for various beam waists in the Maxwell equations 
and become smaller for larger beam waists in the momen- 
tum conservation law. Moreover, refining the discretiza- 
tion helps to reduce both residuals making the dipole- 
based beam model numerically superior with respect to 
the Gaussian approximation when calculating the force 
densities, especially for tightly focused beams. 

Our numerical experiments showed that dipole-based 
beams very closely resemble ideal Gaussian beams for 
large waists and contain expected diffraction artifacts 
(side-lobes) with smaller waists. For example, Figure [2] 
shows the result for two linearly polarized beams with 
different wavelengths and the same desired waist. The 
computational domain is a cube with 10 /im sides. The 
dipole array is situated 2 /im outside the domain and rep- 
resents a 15A x 15A square aperture uniformly filled with 
32 tuned point dipoles in each direction. It is important 
to realize that the dipole model is not always able to 
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FIG. 2: Intensity profiles of two focused beams with the same 
desired waists but different wavelengths modeled via the su- 
perposition of dipoles tuned to mimic ideal Gaussian beams 
(top: Ao = 1.2 pm, bottom: Ao = 3.2 /im). Actual waists are 
shown as contour lines. The larger wavelength beam (bot- 
tom row) corresponds to a smaller relative desired waste (i.e. 
larger beam parameter \s\) and deviates from the Gaussian 
beam while conserving the momentum. Both beams have the 
same power, and the drop in the intensity at Ao = 3.2 /im is 
due to the increased absorption in water. 

achieve the waist size of the original Gaussian beam it 
mimics (see the contour lines in Figure [2J . This is the 
price one pays for a better conservation of momentum. 

The power emitted by the dipole array into the liquid is 
determined and adjusted by computing the integral of the 
normal component of the Poynting vector over the (vir- 
tual) interface between the computational domain and 
the array. Also, in this paper we neglect the reflection of 
light at the interface between the vessel and the surround- 
ing medium, considering the whole medium to be opti- 
cally homogeneous, although the liquid is confined to a 
finite spatial domain. The embarrassingly parallel nature 
of the superposition principle allowed us to compute the 
electromagnetic field and the corresponding force density 
very efficiently by exploiting the computational power of 
the Graphics Processing Unit (GPU). 



IV. LIGHT-DRIVEN INCOMPRESSIBLE FLUID 

The Navier-Stokes equations for an incompressible 
Newtonian fluid are (see e.g. 

V-v = 0, (17) 
Poftv + p (v V)v-//V 2 v + Vp = f, (18) 

where v(x, t) is the local velocity of fluid, p(x, t) is the 
pressure, po is the mass density, p is the viscosity and 
f (x, t) is the applied volumetric body force density. Since 
light is partially absorbed by the liquid, we need to con- 



sider the following heat equation: 

p c p d t T + p c p {v ■ V)T-k\7 2 T = q, (19) 

where T(x, t) is the temperature, c p is the specific heat 
capacity at constant pressure, and q(x, t) is the heat 
source density. 

To model the motion of fluid due to the heat-driven 
expansion we use the Boussinesq approximation, which 
amounts to splitting the body force into two parts: 

f=(f L )+pg, (20) 

where g is the acceleration due to gravity, and (f L ) is the 
time-averaged Lorentz force density ([T31) . The modified 
mass density p(x, t) is described by 

p = po [1 - /3(T - To)} , (21) 

where To is the reference temperature (before the heat 
source q is applied), and f3 is the thermal expansion co- 
efficient of the fluid. 

In the present paper we consider the stationary flows 
only, so that all quantities above are considered to be 
time- independent. This also simplifies the equations (|18[) 
and (1191) . We impose the no-slip boundary condition v = 
at the walls of the vessel. 

The numerical solution of the above coupled system 
of non-linear equations was implemented in the open- 
source software environment OpenFOAM [ll|, which is 
based on the finite- volume discretization scheme and fea- 
tures a rich set of robust algorithms. In particular, we 
employed the iterative SIMPLE (Semi-Implicit Method 
for Pressure-Linked Equations) algorithm [l6|, Oj} that 
converged reasonably fast to the expected tolerance in 
all our numerical experiments. We have also tested the 
convergence of the numerical solution (to a stable result) 
as a function of the discretization step and observed the 
expected second-order behavior. 

V. ANALYSIS OF FLOWS 

In this section we describe the flows that would be in- 
duced in pure water at the reference temperature To = 
300 K if the Lorentz force model was indeed true. We 
have computed the induced flows at several wavelengths 
in the neighborhood of the pronounced absorption peak 
at 3 /an shown in Figure[3l |18j . The wavelengths 1.2 /im, 
2.4 pm, 2.8 pm, and 3.2 /im indicated in Figure [3] cover 
the regions of transparency, as well as normal and anoma- 
lous dispersion allowing to distinguish between the effects 
of the different terms in the Lorentz force density (p~3|) 
and the heat source (|T3|) . The most interesting results 
were obtained at 1.2 /mi where the water is almost loss- 
less and the body force (|13[) is dominated by the first 
term while the heat source (fT5j) is almost zero, and at 
3.2 pm where, although the absorption is significant, it 
does not extinct the beam too soon and lets it propagate 
some length. 



5 





100 fim/s — 




FIG. 3: Real and imaginary parts of the complex relative 
permittivity of water in the neighborhood of the absorption 
peak at 3/^m. Circles indicate different wavelengths (1.2 /im, 
2.4 /im, 2.8 /im, and 3.2 /im) considered in our numerical ex- 
periments. 



The container and the computational domain of our 
simulations is a cube with 10 /um sides. The Gaussian 
beam (re-modeled by the dipole array) propagates along 
one of the Cartesian axes coinciding with the edge of 
the domain and is linearly polarized along one of the 
other edges. We have simulated a whole range of incident 
powers with the results presented below corresponding to 
lmW. 

In all our experiments the spatial pattern of the flow 
appears to be divided into four quarters by the two planes 
intersecting at the beam axis. One of the planes is the 
polarization plane of the beam and the other plane is 
orthogonal to it. Despite the presence of gravity, which 
in our case acts along the polarization direction, i.e., or- 
thogonal to the beam axis, the flow patterns are largely 
symmetric about the mentioned two planes. 

Due to the assumed incompressibility of water the re- 
sulting stationary flows are divergence-free. Since the 
water container is closed and there are no sources or sinks 
the stream lines should be closed. This is exactly what 
we observe in Figure 0] where the stream lines and the 
quiver plot of the induced flow at 1.2 /im are shown. The 
axis of the beam is vertical and its direction of prop- 
agation is upwards in this and subsequent 3D plots. In 
each of the aforementioned quarters the stream lines rep- 
resent simple closed concentric loops centered around a 
line parallel to the beam axis. The direction of flow along 
each such loop is always away from the beam axis along 
the polarization direction and coming back towards the 
beam along the orthogonal direction, see Figure [4] (bot- 
tom) . The velocity along the loops is highest at the point 
closest to the beam axis reaching ~ 87/im/s. 

The flow pattern changes when we consider the wave- 
lengths where the optical absorption in water becomes 
significant introducing the second term in the body force 
(fTB")) . For example, in Figure [S] we show some of the 



FIG. 4: Induced flow at 1.2 /im (low optical losses). Top: typi- 
cal stream lines (beam propagates vertically upwards through 
the middle of the domain). Bottom: quiver plot in the plane 
orthogonal to the beam axis showing the four characteristic 
quarters. 

stream lines of the flow induced at 3.2 fan. There are just 
two simple closed loops in each quarter in this case (both 
are shown). The rest of the stream lines have a more 
complicated shape featuring a large number of windings. 
It is also apparent that the stream lines are no longer 
orthogonal to the beam axis. The direction of flow is de- 
termined by both the polarization and the propagation 
direction of the beam in this case. Although the heat 
term (| 1 5|) is no longer zero and the direction of grav- 
ity breaks the symmetry of the problem, the overall flow 
pattern remains almost symmetric about the same two 
orthogonal planes as in the lossless case. 

Our simulations show that the flows induced by the 
beams with higher powers have the same spatial pattern 
as low-power flows and the computed velocity scales lin- 
early with the body force up to 50 mW. This is a clear 
indication of the low Reynolds number regime, meaning 
that in the future studies a linear Stokes approximation 
can be used to simplify the computational model. 



VI. CONCLUSIONS 

We have demonstrated that the induced Lorentz force 
model (|11[) ~ (I13I) results in significant hydrodynamic flows 
in the neighborhood of a single-beam optical trap in a 
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FIG. 5: Stream lines of the induced flow at 3.2 jj,m (high optical losses). Each plot shows a single stream line seeded at a 
particular location and propagated for the same large number of steps with the Runge-Kutta algorithm. In each quarter only 
two stream lines are simple closed loops (both shown). Other stream lines have a more complex shape (the typical ones are 
shown). The flow patterns in the other three quarters of the domain are similar. 



liquid environment. The existence of such flows would 
suggest the form of the electromagnetic momentum con- 
servation law (J5J) different from the ones due to Abraham 
and Minkowski, which after time averaging both feature 
the Helmholtz force density that equals zero in a homo- 
geneous liquid. Of course, the fact that such flows have 
not been directly observed so far (except, perhaps, in 
[E BE3) could simply indicate their absence. These in- 
duced flows, however, are almost invisible by their nature 
and may be difficult to detect with other means. 

Normally a flow pattern can be visualized using small 
tracer particles, such as ink. This is not straightforward 
in the present case, since both the flow and the tracer 
particle will be influenced by the laser beam. The equa- 
tion of motion for a small spherical particle with mass m 
and position x(t) is 



d 2 x 



(22) 



where F om and Fdrag are the electromagnetic force and 
the fluid drag force acting on the particle. 

Considering a sufficiently small tracer particle without 
optical losses the electromagnetic force can be approxi- 
mated by the gradient force: 



F grad = 7rr 3 Re( ew ) £p / g ° \ ^ 1 

£p/£0 + 2 



V|E| 2 , 



(23) 



where e w and e p are the permittivities of the water 
and particle respectively. The Stokes drag is given by 
Fdrag = 67r/^ w ru, where /i w is the dynamic viscosity of 
water, r is the particle radius, and u is the particle ve- 
locity with respect to the stationary fluid. In the present 
case, however, the fluid is already in motion due to the 



induced flow. Let the velocity of this flow with respect 
to the stationary reference frame be v. Then, the drag 
force is Fd rag = 67r/Lt w r(v — dx/dt), and the equation of 
motion becomes 



d 2 x dx 
m _ + 67rMwr _ 



7rr 3 Re(e w ) 



gp/gp - 1 

e p /e + 2 



(24) 



V|E| 2 + 67^ rv. 



Paradoxically, if the tracer velocity dx/dt coincides with 
the velocity v of the fluid flow, i.e., the tracer serves its 
purpose and moves with the flow, then the drag force dis- 
appears, the tracer will be pulled by the optical gradient 
force and thus no longer follow the flow. On the other 
hand, if the gradient force pulls the particle in the direc- 
tion opposing the flow, the friction term will slow it down 
and the drag force will reduce the trapping efficiency. 

Fortunately, the difference in the dependence on the 
particle radius (r 3 versus r) means that sufficiently small 
tracers will mostly follow the induced flow, since the gra- 
dient optical force acting on them will be much smaller 
than the drag force. However, such tracers will have to 
be much smaller than the wavelength making it difficult 
to actually see their motion. 

Finally, the equation of motion (|24p does not include 
the stochastic Brownian term which will make tracers 
jump from one streamline to another, further complicat- 
ing the detection of the simulated flow patterns. Hence, 
a detailed analysis of the particle dynamics and a ro- 
bust experimental technique are needed before the in- 
duced flows predicted in this paper can be confirmed or 
truly ruled out. 



7 



[1] A. Ashkin, J. M. Dziedzic, J. E. Bjorkholm, and S. Chu, 
Observation of a single-beam gradient force optical trap 
for dielectric particles, Opt. Lett, 11 288-290, 1986. 

[2] M. Dienerowitz, M. Mazilu, and K. Dholakia, Opti- 
cal manipulation of nanoparticles: a review, Journal of 
Nanophotonics, 2, 021875, 2008. 

[3] J. E. Molloy and M. J. Padgett, Lights, action: optical 
tweezers, Contemporary Physics, 43, 241-258, 2002. 

[4] E. J. G. Peterman, F. Gittes, and C. F. Schmidt, Laser- 
induced heating in optical traps, Biophysical Journal, 84, 
1308-1316, 2003. 

[5] R. T. Schermer, C. C. Olson, J. P. Coleman, and F. Bu- 
choltz, Laser-induced thermophoresis of individual parti- 
cles in a viscous liquid, Optics Express, 19, 10571, 2011. 

[6] F. M. Weinert and D. Braun, Optically driven fluid flow 
along arbitrary microscale patterns using thermoviscous 
expansion, J. Appl. Phys., 104, 104701, 2008. 

[7] P. Wu, R. Huang, C. Tischer, A. Jonas, and E.-L. Florin, 
Direct measurement of the nonconservative force field 
generated by optical tweezers, Phys. Rev. Lett., 103, 
108101, 2009. 

[8] M. Khan and A. K. Sood, Tunable Brownian vortex at 
the interface, Phys. Rev. E, 83, 041408, 2011. 

[9] D. Psaltis, S. R. Quake, and C. Yang, Developing optoflu- 
idic technology through the fusion of microfluidics and 
optics, Nature, 442, 381-386, 2006. 



[10] O. A. Louchev, S. Juodkazis, N. Murazawa, S. Wada, 
H. Misawa, Coupled laser molecular trapping, cluster as- 
sembly, and deposition fed by laser-induced Marangoni 
convection, Optics Express, 16, 5673, 2008. 

[11] A. H. J. Yang, T. Lerdsuchatawanich, and D. Erickson, 
Forces and transport velocities for a particle in a slot 
waveguide, Nana Lett, 9(3), 1182-1188, 2009. 

[12] J. C. Ryu, H. J. Park, J. K. Park, and K. H. Kang, New 
electrohydrodynamic flow caused by the Onsager effect, 
Phys. Rev. Lett, 104, 104502, 2010. 

[13] W. Frias and A. I. Smolyakov, Electromagnetic forces 
and internal stresses in dielectric media, Phys. Rev. E, 
85, 046606, 2012. 

[14] G. K. Batchelor, An introduction to fluid dynamics, Cam- 
bridge University Press, Cambridge, 2000. 

[15] OpenFOAM 2.1.1 User Guide, Open CFD Limited, 
http: //www . openfoam.org/docs/user/ , 2012. 

[16] H. Jasak, Error analysis and estimation for the Finite 
Volume method with applications to fluid flows, PhD The- 
sis, Imperial College, University of London, 1996. 

[17] S. V. Patankar, Numerical heat transfer and fluid flow, 
Hemisphere Publishing Corporation, 1980. 

[18] G. M. Hale and M. R. Querry, Optical Constants of Wa- 
ter in the 200-nm to 200-/^m Wavelength Region, Appl. 
Opt, 12(3), 555-563, 1973. 



